Iterative Method#
강좌: 수치해석
Sparse Matrix#
많은 element가 0이고 일부만 값을 가지는 매우 크고 성긴 행렬임
공학 문제에서 매우 큰 Sparse Matrix 해석이 필요함
Sparse matrix를 효과적으로 저장하고, 계산하는 전용 라이브러리가 존재함
scipy.sparse내 다양한 함수가 존재함
Iterative Method#
개념#
매우 큰 행렬 System \(Ax=b\) 를 반복해서 푸는 방법이다.
기본 개념은 다음과 같다.
\(A = A_1 - A_2\)
\(A_1\) 은 역행렬을 쉽게 구해지는 형태이다.
\[ A_1 x = A_2 x + b \]반복되는 해를 \(x^{(k)}\) 하고 이를 적용한다.
\[ A_1 x^{(k+1)} = A_2 x^{(k)}+ b \]\(x^{(k)} \rightarrow x\) 이면 오차 \(e^{(k)} = x^{(k+1)} - x^{(k)} \rightarrow 0\) 이다. 즉 오차가 \(e^{(k)}\) 감소할 때 까지 반복한다.
모든 경우에 오차가 감소하지 않는다. \(A_1^{-1} A_2\) 의 Eigenvalue가 모두 1 보다 작아야 한다.
예제#
지난 시간에 다룬 삼중 대각 행렬을 예제로 한다.
많은 공학 시뮬레이션에서는 Banded Matrix가 많이 사용된다.
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.dpi'] = 150
# Dimension
n = 10
# Matrix system
A = np.diag([2.04]*(n)) + np.diag([-1]*(n-1), -1) + np.diag([-1]*(n-1), 1)
# Forcing term
b = np.ones(n)*0.8
b[0] += 40
b[-1] += 200
# Solve
np.linalg.solve(A, b)
array([ 44.42707508, 49.83123316, 56.42864057, 64.48319361,
74.31707438, 86.32363813, 100.98314741, 118.88198259,
140.73609706, 167.41965542])
Point Jacobi Method#
이 방법은 \(A_1 = D\) 로 한 경우이다.
위 문제의 각 행은 다음과 같다.
Jacobi Method를 적용하면 (n+1) step 에서 i 번째 항은 (n) 번째 Off-diagnal 항을 이용해서 계산할 수 있다.
구현하는 과정은 다음과 같다.
모든 i 에 대해서 \(\Delta T^{(n)}_i\) 를 구함
모든 i 에 대해서 \(T^{(n+1)}_i\) 최신화 함
def jacobi(n, t, b, dt):
"""
Jacobi method
Parameters
----------
n : integer
size
t : array
current solution
b : array
forcing term
dt : array
difference
"""
for i in range(1, n+1):
dt[i] = (t[i-1] + t[i+1] + b[i-1])/2.04 - t[i]
n = 10
tol = 1e-8
# 양 끝점을 포함해서 행렬 만듬
t = np.zeros(n+2)
dt = np.empty_like(t)
# Forcing term
b = np.ones(n)*0.8
b[0] += 40
b[-1] += 200
err = 1
hist_jacobi = []
while err > tol:
# Run Jacobi
jacobi(n, t, b, dt)
# Compute Error
err = np.linalg.norm(dt) / n
hist_jacobi.append(err)
# Update solution
t += dt
# Solution
print(t[1:-1])
[ 44.42707492 49.83123285 56.42864015 64.48319308 74.31707383
86.32363756 100.9831469 118.88198215 140.73609676 167.41965526]
# Validation
A @ t[1:-1]
array([ 40.79999999, 0.79999994, 0.79999998, 0.7999999 ,
0.79999997, 0.79999989, 0.79999998, 0.79999992,
0.79999999, 200.79999997])
plt.semilogy(hist_jacobi)
plt.legend(['Jacobi'])
plt.xlabel('Iter')
plt.ylabel('Error')
Text(0, 0.5, 'Error')
Gauss-Seidel#
Jacobi Method를 생각하면, (n+1) step에서 i번째 값을 구할 때 i-1번째 값의 (n+1) 번째 값을 알 수 있다.
즉 위 식을 다음과 같이 생각할 수 있다.
행렬 식으로 생각하면 \(A=L+D+U\) 라 생각했을 때 \(A_1 = L + D\)인 방법이다.
구현하는 과정은 loop를 하나로 합친다.
i 번째 값에 대해 \(T^{(n+1)}_i\) 최신화 함
i 번째 값에 대해 \(\Delta T^{(n)}_i\) 를 구함
def gauss_seidel(n, t, b, dt):
"""
Jacobi method
Parameters
----------
n : integer
size
t : array
current solution
b : array
forcing term
dt : array
difference
"""
for i in range(1, n+1):
ti = t[i]
t[i] = (t[i-1] + t[i+1] + b[i-1])/2.04
dt[i] = t[i] - ti
n = 10
tol = 1e-8
# 양 끝점을 포함해서 행렬 만듬
t = np.zeros(n+2)
dt = np.empty_like(t)
# Forcing term
b = np.ones(n)*0.8
b[0] += 40
b[-1] += 200
err = 1
hist_gs = []
while err > tol:
# Run Jacobi
gauss_seidel(n, t, b, dt)
# Compute Error
err = np.linalg.norm(dt) / n
hist_gs.append(err)
# Solution
print(t[1:-1])
[ 44.42707497 49.83123296 56.42864031 64.48319331 74.31707408
86.32363785 100.98314716 118.88198239 140.73609693 167.41965536]
plt.semilogy(hist_jacobi)
plt.semilogy(hist_gs)
plt.legend(['Jacobi', 'Gauss-Seidel'])
plt.xlabel('Iter')
plt.ylabel('Error')
Text(0, 0.5, 'Error')
비고#
Iterative 기법은 매 반복 횟수마다 \(A_1^{-1} A_2\) 행렬의 Eigenvalue 만큼 오차가 감쇄함
\[\begin{split} \begin{align} x^{(k+1)} &= (A_1^{-1}A_2) x^{(k)} + A_1^{-1} b \\ x^{(k)} &= (A_1^{-1}A_2) x^{(k-1)} + A_1^{-1} b \\ \Delta x^{(k+1)} &= (A_1^{-1}A_2) \Delta x^{(k)} \end{align} \end{split}\]Jacobi와 Gauss Seidel 기법 Eigenvalue 비교
# Jacobi
A1 = np.diag([2.04]*(n))
A2 = A - A1
# Eigenvalue, Eigenvector
eig = np.linalg.eig(np.linalg.inv(A1) @ A2)
l_jac = sorted(np.abs(eig[0]))
print(l_jac)
[0.1395243512483186, 0.13952435124831886, 0.40726962059008465, 0.4072696205900847, 0.6420203273973373, 0.6420203273973383, 0.8247583655207655, 0.8247583655207658, 0.9406793858965665, 0.9406793858965665]
# GS
A1 = np.diag([2.04]*(n)) + np.diag([-1]*(n-1), -1)
A2 = np.diag([-1]*(n-1), 1)
# Eigenvalue, Eigenvector
eig = np.linalg.eig(np.linalg.inv(A1) @ A2)
l_gs = sorted(np.abs(eig[0]))
print(l_gs)
[0.0, 7.154396559720807e-05, 7.16105441542434e-05, 7.16105441542434e-05, 7.167790710375776e-05, 0.019467044587691305, 0.16586854385559155, 0.41219010079138535, 0.6802263614964851, 0.8848777070507405]
# Jacobi 10회 반복, GS 5회 반복
l_jac[-1]**10, l_gs[-1]**5
(0.5425206454135214, 0.5425206454135184)
Succesive overelaxation (SOR)#
Gauss-Seidel 결과에 relxation parameter \(\omega\) 를 이용해서 계산 결과를 가속 또는 감속할 수 있다.
즉
여기서 \(\tilde{T}_{i}^{(n+1)}\) 은 Gauss-Seidel 결과이다.
\(\omega \in 0, 2)\) 에 따라 감속 또는 가속된다.
\(0 < \omega < 1\) - Gauss-Seidel에 비해 감속
\(\omega = 1\) - Gauss-Seidel 방법
\(1 < \omega < 2\) - Gauss-Seidel에 비해 가속
최적의 값은 문제마다 달라지므로 찾아야 한다.
DIY#
위 예제에 대해서 i번째 값에 대한 SOR 수식을 유도하시오.
\[ T^{(n+1)}_i = \alpha_{i-1} T^{(n)}_{i-1} + \alpha_{i} T^{(n)}_{i} + \alpha_{i+1} T^{(n)}_{i+1} + \beta b_i \]SOR 방법을 구현하시오. \(\omega \in (1.1, 1.8)\) 사이에서 값을 변화시키면서 수렴 속도를 비교하시오.
Point Jacobi, Gauss Seidel, SOR 방법에 대해 격자 크기를 달리하면서 해의 변화를 관찰하고, 계산시간의 증가를 파악하라.